Constraints on f{R) gravity from probing the large-scale structure 
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We study cosmological constraints on metric f{R) gravity models that are designed to reproduce 
the ACDM expansion history with modifications to gravity described by a supplementary cosmo- 
logical freedom, the Compton wavelength parameter Bq. We conduct a Markov chain Monte Carlo 
analysis on the parameter space, utilizing the geometrical constraints from supernovae distances, 
the baryon acoustic oscillation distances, and the Hubble constant, along with all of the cosmic 
microwave background data, including the largest scales, its correlation with galaxies, and a probe 
of the relation between weak gravitational lensing and galaxy flows. The strongest constraints, 
however, are obtained through the inclusion of data from cluster abundance. Using all of the data, 
we infer a bound of Bq < 1.1 x 10"^ at the 95% C.L. 



I. INTRODUCTION 

Cosmic acceleration can either be explained by intro- 
ducing large amounts of dark energy or considering mod- 
ifications to gravity such as the addition of a suitable 
function f(R) of the Ricci scalar to the Einstcin-Hilbcrt 
action [iHSl- In fact, one may interpret the cosmologi- 
cal constant as being of this kind rather than attributing 
it to vacuum energy. It has been argued that a valid 
/(i?) model should closely match the ACDM expansion 
history We specialize our considerations to func- 

tions f[R) that exactly reproduce this background and 
parametrize the class of solutions in terms of its Compton 
wavelength parameter Bq Q. Such /(i?) modifications 
affect gravity at solar-system scales, which are well tested 
and impose stringent constraints on deviations from gen- 
eral relativity. However, the chameleon effect 0-13 pro- 
vides a mechanism that allows certain f{R) models to 
evade solar-system tests (e.g., Q)- The transition re- 
quired to interpolate between the low curvature of the 
large-scale structure and the high curvature of the galac- 
tic halo sets the strongest bound on the cosmological field 

(So < 10-5 a El). 

Independently, strong constraints can also be inferred 
from the large-scale structure alone. The enhanced 
growth of structure observed in f[R) gravity mod- 
els manifests itself on the largest scales of the cosmic 
microwave background (CMB) temperature anisotropy 
power spectrum @, where consistency with CMB data 
places an upper bound on Bq of order unity Cross 
correlations of the CMB temperature field with fore- 
ground galaxies serve as another interesting test of f{R) 
gravity models [ll| , tightening the constraint on the 
Compton wavelength parameter by an order of mag- 
nitude (e.g., [13 )• However, the currently strongest 
constraints on f{R) gravity from large scale structures 



are inferred from the analysis of the abundance of low- 
redshift X-ray clusters, yielding an improvement over the 
CMB constraints on the free field amplitude of the Hu- 
Sawicki l4| [n = 1) model of nearly four orders of mag- 
nitude [13 1 . 

In this paper, we perform an independent analysis of 
f{R) gravity constraints, focusing on cosmological data 
only. Our analysis differs from previous studies par- 
tially in terms of the theoretical model, the parametrical 
approach, or data sets implemented. We compare and 
discuss our results to the constraints of previous anal- 
yses. Here, we conduct a Markov chain Monte Carlo 
(MCMC) study of metric f(R) gravity models that are 
designed to reproduce the ACDM expansion history us- 
ing data from CMB anisotropics, supernovae distances, 
the baryon acoustic oscillation (BAO) distances, and the 
Hubble constant. For observables in the linear regime, we 
adopt the parametrized post-Friedmann (PPF) frame- 
work [l5| and its implementation into a standard 
Einstein-Boltzmann linear theory solver [ll for the the- 
oretical predictions. This framework allows us to include 
information from the near horizon scales. We also uti- 
lize information from the cross correlation between high- 
redshift galaxies and the CMB through the integrated 
Sachs- Wolfe (ISW) effect. We further use a probe of the 
relation between weak gravitational lensing and galaxy 
flows, as well as data from the abundance of clusters 
that are identified by overdensities of bright, uniformly 
red galaxies. Latter yields the tightest constraints on the 
cosmological parameters, particularly on Bq. We com- 
pare these constraints to the results of [H derived for 
the Hu-Sawicki model. 

In we review metric /(i?) gravity theory. We 
present the results of our MCMC study in mil] and dis- 
cuss them in ijlVI Finally, the PPF formalism for f{R) 
gravity [31 and details about the modifications to the 
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ISWWLL code [13, li3 used for the galaxy-ISW (gISW) 
cross correlation observations are specified in the Ap- 
pendix. 



II. f[R) GRAVITY 

In f{R) gravity, the Einstein-Hilbert action is supple- 
mented by a free function of the Ricci scalar R, 



This can be solved numerically with the initial conditions 



(1) 



Here, = SttG and £m is the matter Lagrangian, where 
we have set c = 1. Variation with respect to the metric 
gap yields the modified Einstein equation for metric f{R) 
gravity, 

G^.+fRR^,- (^^ - Ufr^ g^.-V,,V,/H = k'T^,, (2) 

where subscripts of R indicate differentiation with re- 
spect to the Ricci scalar and the connection is of Levi- 
Civita type. The modified Friedmann equation is derived 
in the usual way, i.e., taking the time-time component of 
Eq. ©, 



- fniHH' + H'-) 



^+H\fnnR' ^^p, (3) 



where here and throughout the paper primes denote 
derivatives with respect to In a. 



/(In ai) 
/'(In a,) 



(6) 
(7) 



where p = (-7 + V73)/4 and 0.01 [g. A is an 

initial growing mode amplitude and characterizes a spe- 
cific solution in the set of functions /(i?) that recovers 
the ACDM background. Note that the amplitude of the 
decaying mode is set to zero in order to not violate con- 
straints in the high-curvature regime. We follow and 
parametrize our solutions in terms of the Compton wave- 
length parameter 



B 



l + fR H' 



(8) 



evaluated at Bq = B{\na = 0) rather than by A. In a 
Taylor expansion of /(i?), the first term corresponds to a 
cosmological constant, while the second is a rescaling of 
Newton's constant. With its fan term, this parametriza- 
tion therefore captures the essence of the modification. 
Standard gravity is recovered in the case where Bq = 0. 
For stability reasons the mass squared of the scalar field 
fn must be positive, which implies Bq > [gI, [20j. 



III. COSMOLOGICAL CONSTRAINTS 



1. Designer model 

The function f{R) can be constructed in a way that it 
recovers any background history. Given the equation of 
state of the effective dark energy, wde, one integrates the 
modified Friedmann equation to obtain the correspond- 
ing form of f{R) @, It has been pointed out that 
viable f{R) cosmologies must closely match the ACDM 
expansion history 0, Here, we focus on flat models 
and consider modifications to gravity that reproduce the 
ACDM expansion history exactly, i.e., wde = ~1j where 
we neglect contributions from radiation. Note that the 
restriction to the matter-dominated epoch is well mo- 
tivated by requiring that the well-tested high-curvature 
regime reproduces standard phenomenology (e.g., [B])- 
The ACDM background is then given by 



= ^(Pm + PA). 



(4) 



Equating it with the matter-dominated Friedmann equa- 
tion, Eq. ([3]), yields an inhomogeneous second order dif- 
ferential equation for f{R), 



f" 



H' R" 



/' + 



R' 
6H' 



R' 



(5) 



We use a variety of cosmological data sets to con- 
strain the f{R) theory of gravity. First we use the CMB 
anisotropy data from the five-year Wilkinson Microwave 
Anisotropy Probe (WMAP) [2l[, the Arcminute Cosmol- 
ogy Bolometer Array Receiver (ACBAR) [2^ . the Cos- 
mic Background Imager (CBI) [23], and the Very Small 
Array (VS A) [23| • Next we employ data from the Super- 
nova Cosmology Project (SCP) Union [2^ compilation, 
the measurement of the Hubble constant from the Su- 
pernovae and Hq for the Equation of State (SHOES) J2y| 
program, and the BAO distance measurements of [27|. 
Furthermore, we take gISW cross correlation observa- 
tions using the iswwll code of [13, [ill , the Eq mea- 
surement, probing the relation between weak gravita- 
tional lensing and galaxy fiows, of 28], as well as cluster 
abundance constraints from the likelihood code of [2^ 
(CA). We quote results with the latter three constraints 
separately to highlight their impact on the results. For 
comparison, wc also analyze constraints that include the 
(Ts measure of the Chandra Cluster Cosmology Project 

(CCCP) [13. 

In mil Al we discuss the predictions for some of these 
observables for specific values of the Compton wavelength 
parameter. In ^III Bl wc present the results of a MCMC 
likelihood analysis, which is conducted with the publicly 
available COSMOMC [111 package. 
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A. Cosmological observables 

In this section we illustrate model predictions of the 
various cosmological observables we use in the con- 
straints. We choose the parameters of the various models 
that highlight results from the MCMC analysis. 

By construction, at high redshifts, the ,f{R) modifica- 
tions become negligible, and so we choose a parametriza- 
tion that separates high-redshift and low-redshift con- 
straints. Specifically we take 6 high-redshift parameters: 
the physical baryon and cold dark matter energy density 
flbh^ and flch^, the ratio of sound horizon to angular di- 
ameter distance at recombination multiplied by a factor 
of 100 9, the optical depth to reionization t, the scalar 
tilt Us, and amplitude As at fc* = 0.002 Mpc~^. 

Since restricting to flat universes, at low redshifts, 
ACDM has no additional parameter, whereas in f{R) 
gravity, an extra degree of freedom is introduced by the 
Compton wavelength parameter Bq , where ACDM is re- 
produced in the limit Bq — > 0. 

We illustrate predictions from the maximum likelihood 
ACDM model using aU of the data (see Table H]). For 
f{R) gravity, where not otherwise specified, we use these 
ACDM best-fit parameters and add nonzero values of Bq 
for comparison. 



1. Geometrical measures 

The comparison of the magnitudes of high-redshift to 
low-redshift supernovae yields a relative distance mea- 
sure. The acoustic peaks in the CMB, the measurement 
of the local Hubble constant, and the BAO distances ad- 
ditionally provide absolute distance probes which com- 
plement the relative distance measure of the supernovae. 
These probes constrain the background and since our 
f{R) models are designed to match the ACDM expan- 
sion history with vanishing effect in the high-curvature 
regime, they do not generate any tension between the 
models. For the Hubble constant, we use the SHOES 
measurement of 



Hq = 74.2 ± 3.6 km s"^ Mpc"\ 



(9) 



which employs Cepheid measurements to link the low- 
redshift supernovae to the distance scale established by 
the maser galaxy NGC 4258. Further, we apply the con- 
straint 



= (0.282 ±0.018) 



0.1326 



0.58 



(10) 



from the BAO distance measure of [27[ that is obtained 
from analyzing the clusterin g of galaxies from the Sloan 
Digital Sky Survey (SDSS) [H and the 2-degree Field 
Galaxy Redshift Survey (2dFGRS) fs^- This measure- 
ment yields the tightest constraint on fl^ and substan- 
tially assists in breaking degeneracies with Bq in the clus- 
ter abundance data. 




a. 4000 



FIG. 1: Best-fit CMB temperature anisotropy power spec- 
trum for ACDM and f{R) gravity from using the Union, 
SHOES, BAO, and CMB data. 



2. The cosmic microwave background 



The CMB probes the geometry of the background his- 
tory as well as the formation of large-scale structure. The 
latter manifests itself on the largest scales through the 
ISW effect from the evolution of the gravitational po- 
tential. To predict these effects we implement the PPF 
modifications described in Appendix [X] The incorpo- 
ration of the PPF formalism into a standard Einstein- 
Boltzmann linear theory solver yields an efficient way to 
obtain predictions of f{R) gravity for the CMB. We uti- 
lize the PPF modifications to CAMB [U implemented 
in Ref. [lB|jWhich we configure for f{R) gravity as de- 
scribed by |15| . 

In Fig. [1] we plot the CMB temperature anisotropy 
power spectrum with respect to angular multipole £ for 
the best-fit ACDM model and the best-fit f {R) grav- 
ity model using the Union, SHOES, BAO, and CMB 
data jointly (see Tables [11 and IIIip . Adding an additional 
freedom from /(i?) gravity only yields an insignificant 
improvement in the fit. Relative to ACDM, or equiv- 
alcntly Bq = 0, the growth of structure is enhanced. 
The ISW effect at the lowest multipoles is decreasing 
for Bo < 3/2 [1, [Il|. At - 3/2, the ISW contribu- 
tion rises again and turns into a relative enhancement 
over ACDM for Bq > 3, ruling out So > 5 by WMAP 
data 

Due to suppression of f{R) modifications in the high- 
curvature regime, the CMB acoustic peaks can be uti- 
lized as usual to infer constraints on the high-redshift 
parameters, in particular the physical energy densities of 
baryonic matter and dark matter as well as the angular 
diameter distance to recombination. 
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3. Galaxy-ISW cross correlations 

The correlation between galaxy number densities and 
the CMB anisotropics can be used to isolate the ISW ef- 
fect in the CMB. Whereas the ISW effect in the CMB is 
suppressed for Bq <3 and enhanced for Bq >3 with re- 
spect to Bq = 0, the gISW cross correlation is suppressed 
for all i?o > and leads eventually to anti-correlations 
(see Fig. [5]). This makes the gISW cross correlation in- 
teresting for improving constraints (see, e.g., @, [HI). As 
was shown in the absence of negative correlations 
between the CMB and an assortment of galaxy surveys 
infer a boundary of i?o ^ 1- This constraint can be im- 
proved by a more rigorous analysis of gISW cross cor- 
relations as was performed by, e.g., This study 
implements a parametrization for the modifications in- 
duced by /(i?) gravity that is based on the introduction 
of an effective scalar degree of freedom in the Einstein- 
Hilbcrt action [H, [s^ and uses gISW cross correlation 
data from (37| . Here, we conduct an independent anal- 
ysis, utilizing the PPF framework and for the calcula- 
tion of the gISW likchhood, the iswwll code of [l7l.[l8|. 
which has turned out to be very useful for constraining 
infrared modifications of gravity j38l - l4]| . 

We evaluate the gISW cross correlation in the Lim- 
ber and quasistatic approximation, as it is done in the 
ISWWLL code used for the data analysis. The gISW cross 
correlation in this approximation reads 



C: 



a + 1/2)2 



dzfj{z)H{z) 



D{a,k) — G{a,k)P{k) 



(11) 



with z = 1/a — 1 denoting redshift. Here, D{a,k) is 
the linear growth rate in the quasistatic regime defined 
by Arn{a, k) = A„,(l, k)D{a, k)/D{l, k), where A„i(a, k) 
is the matter density perturbation. P{k) is the matter 
power spectrum today. 

The Limber approximation becomes accurate at the 
percent level for £ > 10. This condition is satisfied by 
about 90% of the total 42 data points that are used 
in the iswwll code. The data are divided into nine 
galaxy sample bins j, i.e., 2MASS0-3, LRGO-1, QSOO-1, 
and NVSS. The function fj{z) relates the matter den- 
sity to the observed projected galaxy overdensity with 
fj{z) = bj{z)Ilj{z) in the absence of magnification bias. 
Hj(z) is the redshift distribution of the galaxies and the 
bias factor bj{z) is assumed independent of scale (cf. j4^). 
but dependent on redshift. The code determines fj{z), 
among other things, from fitting auto power spectra and 
cross power spectra between the samples. We refer to 
Appendix |B] for more details about the data and the ac- 
curacy of the Limber approximation. 

Scalar linear perturbations of the Friedmann metric 
are presented here in longitudinal gauge, i.e.. 



where dx"^ is the unperturbed spatial line element with 
curvature K = 0. We define $_ = ($- *)/2. In the 
quasistatic regime, we infer 



G{a,k) = 



$-(a,fc) A^{ai,k) 1 
$_(«;, fc) Am(l,fc) a^ 
1 D{a,k)l 



l + fRD{l,k)a 
from the modified Poisson equation. 



(13) 



2 1 + //? a 



A„,(a,fc), (14) 



requiring that the initial conditions lie well within the 
high-curvature regime, where general relativity holds. 
This implies <i>_(aj;,fc) ~ $_(ai), D{ai,k) ~ D{ai), and 
fn^ai) ~ 0. We solve the ordinary differential equation 



A:: + (2 + aA'_^i^^SVA. 



H 



2 1 + fR H^a^ 



(15) 



for the linear matter density perturbation Am(a, fc), 
where g{a,k) is the metric ratio in PPF formalism (see 
Appendix Note that in the limit Bq — > 0, we have 
g — > 0. Therefore, in this limit, Eq. ([T5|) recovers the 
quasistatic ordinary differential equation for the matter 
overdensity in ACDM. Wc solve Eq. (|15p with initial con- 
ditions at Oi 1, in a regime where general relativity 
is expected to hold, i.e., A'^{ai,k) — A,„(ai,fc) with a 
normalization set by the initial power spectrum. At the 
scales that are relevant for the gISW cross correlations, 
the product D dG/dz used in Eq. ([TT|) is accurately de- 
scribed through solving Eq. ([T5|) for Am and using the 
approximated G{a, k) of Eq. p^ . We show this by com- 
paring the approximated product to its counterpart from 
a full derivation within linear perturbation theory (see 
Fig. [3]). We take the relations that exactly describe the 
scalar linear perturbation theory in f{R) gravity for a 
matter-only universe from Ref. 



4-. Weak gravitational tensing and galaxy flows 

The relationship of weak gravitational lensing around 
galaxies to their large-scale velocities has been proposed 
as a smoking gun of gravity (isj . The advantage of such 
a probe lies in its insensitivity to galaxy bias and initial 
matter fluctuations. The expectation value of the ratio 
of galaxy-galaxy to galaxy- velocity cross correlations of 
the same galaxies yields an estimator [i^ 



Er 



(l + ./i?)/3 



(16) 



(1 + 2^)dt^ + a\l + 2^)dx^ 



(12) 



where (3 = dlnAm/rflna. Recently this quantity has 
been measured analyzin g 7 205 luminous red galax- 
4l from the SDSS [H, yielding Eg = 0.392 ± 
by averaging over scales R — (10— 50)/i^^ Mpc. 
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increasing effective, bias-weighted, redshift. Adding nonzero values for the Compton wavelength parameter (red dashed lines): 
Bo = 0.1 (top) and Bo = 0.5 (bottom). 



Note that the Eg measurement is based on spectroscopic 
LRG samples, whereas the gISW analysis uses photomet- 
ric LRG samples that do not overlap in redshift. Further- 
more, the error on Eg is dominated by uncertainties in 
lensing and redshift space distortions and most of the 
signal comes from small scales around lOh^^ Mpc. The 
gISW signal is dominated by large scales and most of 
the error is caused by sampling variance and shot noise 
of galaxies. Therefore, we can safely neglect correlations 
between the Eg and gISW data sets. 

We calculate from solving Eq. ([15]), which yields a 
good approximation to Eg for the scales and Compton 
wavelength parameters of interest. We illustrate predic- 
tions for Eg in Fig. 21 The red-dashed lines show the 
approximated values and crosses indicate check points 
derived using full linear perturbation theory. At the red- 
shift of the measurement, z ~ 0.32, and in the linear 
regime of f{R) gravity, Eq. (|16|) is only weakly dependent 
on scale and shows no fc-dependence at all for Bq = 0. 
Therefore, we only need to evaluate Eg at a representa- 
tive scale, which we choose here as fc = O.lh Mpc~^. We 
can then compare this value to the mean Eg observation 
from [23. At small values of the Compton wavelength 



parameter {Bq < 0.01), a fc-dependence shows up for 
k < O.lh Mpc^^. However, the effect is small compared 
to the error in Eg and is subdominant to the uncertainty 
in Oni. 

Fig. m demonstrates that the Eg probe has the po- 
tential to discriminate between ACDM and f{R) gravity 
even at low values of Bq ■ 



5. Cluster abundance 

Tighter constraints on modified gravity theories and 
in particular on f{R) gravity than from the linear theory 
probes can be achieved by testing the weakly to fully non- 
linear scales (cf. [l^lH, We use the abundance of 
clusters to set our strongest boundary on the possible Bq 
values. We employ a preliminary version of the likelihood 
code of [2^, which utilizes constraints from the most 
massive halos (M > 10^^ Mq//i) inferred from SDSS 
data. The galaxy-galaxy lensing signal from clusters and 
groups from the MaxBCG catalog [i^ is measured in 
three cumulative mass bins corresponding to the nomi- 
nal number densities of groups of 2.5 x 10~^ {Mpc/h)~^, 



6 



0.6 



0.2 - 



0.0 



-0.2 
0.0 



From top: 


So = 0.1 


k = 0.1Ho 




k = 1 //o 




k = 10 Hq 




k = 100 Hq 











0.6 



0.4 



- a 0.2 



0.0 



0.2 



0.4 



0.6 



1.0 



-0.2 
0.0 



hrom top: 


Bo = 0.5 




k = 0.1i/o 






k = 1 //o 






k = 10 Ho 






k = 100 Hq 








^^^^^ 


- 













0.2 



0.4 



0.6 



1.0 



0.6 



0.4 



0.2 - 



0.0 



-0.2 
0.0 



From top: k = 10 //o 






Bo = (ACDM) 






B„ = 0.1 






Bo = 1.0 






Bo = 1.5 













0.6 



0.4 



- S 0.2 



From top: k = 100 //o 

Bo = (ACDM) 
Bo = 0.1 
Bo = 1.0 
Bo = 1.5 



0.2 



0.4 



0.6 



1.0 




FIG. 3: Product of the linear density growth rate D and the derivative of the linear potential growth rate G with respect 
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FIG. 4: Best-fit ACDM (Bo = 0) Ea prediction from linear 
theory (blue solid line). Left panel: adding nonzero values 
for the Compton wavelength parameter (red dashed lines): 
Bo = 1 (top) and Bq = 0.1 (bottom). Right panel: Ea at 
different values of Bo, evaluated at A; = O.l/i Mpc~^. 



2 X (Mpc//i)~3^ and 1.8 x IQ-^ (Mpc//i)-3 and in 

two redshift bins of z = 0.18 and z = 0.25. This signal 
is compared to the theoretical predictions based on the 
mass function calibrated with iV-body simulations, cor- 
rectly taking into account mass-observable scatter, cali- 
bration uncertainties and covariances. We refer to Ap- 
pendix [C] for more details. 

For comparison, we additionally ai ialy ze constraints 
from applying the erg measurement of [30| 



(^8 



0.25 



0.47 



= 0.813 ± 0.013 (stat) ± 0.024 (sys) (17) 



inferred from Chandra observations of X-ray galaxy clus- 
ter samples detected in the ROSAT All-Sky Survey by 
normalizing the mass function at low redshifts. The ha- 
los in the sample have masses AI > 10^^ Mq/H. This 
data was used by jl^ to put a constraint of Ifm] < 
{l.StoD ^ 10"'' (^0 < 10"^) (95% C.L.) on the Hu- 
Sawicki [1] f{R) gravity model. The range indicates a 
=F9% mass calibration error corresponding to the sys- 
tematic error in Eq. ([TT]). Note that in the likelihood 
analysis, we add the systematic error in quadrature. The 
constraint is obtained from using modified forces in the 
spherical collapse calculations. Using standard spherical 
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collapse instead reduces this number to |//jo| < 0.4x 10~^ 
(95% C.L.) [l^. An estimate for the range induced from 
the mass calibration error can be obtained through scal- 
ing from the former result, i.e., Ifml < (0-4lo;2) x 10""'. 
This corresponds to an upper bound of Ifm] < 1.4 x 10~^ 
at the 95% C.L. of both the statistical and systematic er- 
ror. 



B. Constraints 

Our basic cosmological parameter set we use in the 
MCMC analysis is P = {rj^/i^, 6*, r, n,, ln[10i°A,] }. 
On this set we implement the following flat priors: 
ilbh^ e (0.01,0.1), ilch^ e (0.045,0.99), 9 e (0.5,10), 
T e (0.01,0.8), Us £ (0.5,1.5), and ln[10i%] G (2.7,4). 
For f{R) gravity models, where P ^ PU {Bo}, we set a 
flat prior on the free Compton wavelength parameter of 
-Bo € (0,10). 

We begin with the analysis of ACDM in Tables U and 
im We show constraints of using separately gISW, Eq, 
CA, and Eg & CA data together with the CMB, SHOES, 
BAO, and Union measurements, as well as when using all 
of the data jointly. We also quote maximum likelihood 
parameters and value. Horizontal lines divide the chain 
parameters from the derived parameters and the best-fit 
(maximum) likelihood. Notice the improved constraints 
for ACDM that arc obtained by the inclusion of the CA 
data. This analysis sets the baseline to which adding the 
/(i?) degree of freedom should be measured. 

In this paper, we study fiat metric /(i?) gravity models 
that reproduce the ACDM expansion history. Therefore, 
the SHOES, BAO, and Union measurements only fix the 
background and do not distinguish standard from mod- 
ified gravity. However, they contribute to the breaking 
of degeneracies that show up in other data sets and help 
tightening the constraints. When introducing the CMB 
probes, some additional information becomes available. 
In fact, we observe that f{R) gravity yields a slightly 
better fit than flat ACDM, which can be attributed to 
the lowering of the temperature anisotropy power spec- 
trum at small £ (see Fig. [1} , but the improvement in the 
fit is not at a significant level (see Table IIIip . In con- 
trast to fiat ACDM, where the inclusion of the gISW 
data does not yield noticeable improvement on the pa- 
rameter constraints 



17| . in the case of f{R) gravity, an 



order of magnitude improvement on the Bq constraint 
is achieved, i.e., Bq < 0.42 at the 95% C.L. This con- 
straint is in perfect agreement with the independent re- 
sult of [l2j , who found an upper bound on the Compton 
wavelength parameter of Bq < 0.4 at the 95% C.L., us- 
ing gISW cross correlations data from [s^l- Including 
the cluster abundance data instead yields another two 
orders of magnitude improvement over the gISW con- 
straint, i.e., Bq < 3.33 x 10~^. Here and throughout this 
section, we quote constraints at the 95% C.L. If we add 
either the gISW or Eg data, this bound tightens by a 
factor of 1.9 or 2.2, respectively. The joint constraint on 




all data 

with CA 

with gISW 
(Bo ^ lO-^Bo) 



0.000 



0.002 



0.004 0.006 
Bo y 10"- Bq 



0.008 



0.010 



FIG. 5: Marginalized likelihood for Bo when using WMAP5, 
ACBAR, CBI, VSA, Union, SHOES, and BAO in combina- 
tion with the additional data sets. For gISW, the Compton 
wavelength parameter is rescaled as Bo — s- 10~'^i3o, i-e., the 
constraint is a factor of 10^ weaker than illustrated. The hor- 
izontal lines indicate 68%, 95%, and 99% confidence levels. 



the Compton wavelength parameters of Bq < 1.12 x 10^"^ 
(a factor of 3.0) (see Table ITVl is our main result and is 
inferred from combining all of the data sets. We therefore 
found that the gISW and Eg probes arc good comple- 
mentary tests to cluster abundance. In Fig. [5] and Fig. [51 
we plot the marginalized likelihood for Bq and the 2D- 
marginalized contours for Bq and f2m, respectively, with 
different combinations of the data sets, illustrating this 
point. It applies particularly to Eq, where bounds on Bq 
become looser in the absence of data from cluster abun- 
dance. This is what one expects from the trend seen in 
Fig. m i.e., for Bq > 0.1 the Eg prediction approaches 
its ACDM value and eventually overshoots it. 

If we additionally include CCCP to our joint set of 
data, we obtain a constraint of Bq < 0.96 x 10""^ or 
equivalents |/flo| < 1-65 x 10"''. If we take CCCP 
with neither gISW, Eq, nor CA the constraint becomes 
Bq < 1.83 X 10-3 or equivalently IJrq] < 3.12x10"^. The 
constraint is a factor of 2.2 times weaker than what we 
estimated from the result of [Tsj when adding the system- 
atic error from mass calibration in the case of standard 
spherical collapse (see Sec. IIII A~5|) . Note that we have 
used cTg predicted directly by ,f{R) gravity in Eq. pT|) 
rather than constraining the rescaled normalization (T§^ 
in ACDM that matches the halo mass function in f{R) 
gravity at the pivot mass Mes = 3.667 x lO^''' Mq/H 
{M = M500) (cf. [13 )■ In the case of standard spherical 
collapse, this leads to an overestimation of f{R) gravity 
effects on erg of ~ 2% for Bq = 10"^. The error of the 
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0.010 




FIG. 6: Contours of 2D marginalized 68%, 95%, and 99% 
confidence boundaries using WMAP5, ACBAR, CBI, VSA, 
Union, SHOES, and BAO in combination with the additional 
data sets. For gISW, the Compton wavelength is rescaled as 
-Bo — >■ 10~'^i3o, i-e., the constraint is a factor of 10^ weaker 
than illustrated. gISW cross correlations favor higher values 
of Qn-i and break the degeneracy between Q,n-i and Bq seen in 
the CA data. 



CCCP measurement in Eq. ([T7| is ~ ±3%. Also note 
that the CA data infer a preliminary constraint of [2^ 



0-8 



0.25 



0.844 ±0.036, 



(18) 



which is a slightly larger value than in Eq. ([T7)) and there- 
fore admits larger values of Bq, i.e., using CCCP instead 
of CA (without gISW and Eg) tightens the 95% C.L. 
boundary on Bq by a factor of 1.8. It is also important 
to keep in mind that the constraint of [Tsj is inferred for 
the Hu-Sawicki model that exhibits an enhanced linear 
growth with respect to the designer model studied here 
(see, e.g., [H). 

Finally, note that since ACDM is reproduced in the 
limit Bq 0, the slightly poorer fits of f{R) gravity 
with respect to ACDM (see Tables Im] and [W]) have to 
be attributed to sampling errors in the chains. 



utilized all of the CMB data, including the lowest mul- 
tipoles, its correlation with galaxies, the comparison of 
weak gravitational lensing to large-scale velocities, and 
the abundance of clusters. 

We report a constraint on the Compton wavelength pa- 
rameter of Bq < 1.1 X 10^"^ at the 95% C.L. from using all 
of the measurements. This result is substantially driven 
by data from cluster abundance. However as the data 
improve, the limits will saturate due to the chameleon 
effect in massive haloes. gISW measures in combination 
with the CMB, supernovae, BAO distance, and Hubble 
constant probes, yield a constraint of Bq < 0.42 (95% 
C.L.), which is an order of magnitude improvement over 
using the CMB alone as probe of the growth of large-scale 
structure in combination with the geometrical measures. 
This highlights the power of gISW measurements as a 
linear theory probe to constrain infrared modifications 
of gravity. 

The Eg measurement of the relationship between weak 
gravitational lensing and galaxy flows does not improve 
bounds on f{R) gravity on its own. However, when used 
as a complementary probe to cluster abundance, it con- 
tributes substantially to our constraints. This can be at- 
tributed to the slow convergence of its prediction toward 
ACDM when — > 0. It is likely that with improved 
data, the Eg probe will become an important discrimi- 
nator for gravity models. 
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Appendix A: PPF linear theory 



IV. DISCUSSION 

We have performed a MCMC analysis on metric f{R) 
gravity models that exactly reproduce the ACDM ex- 
pansion history. In addition to geometrical probes from 
supernovae, BAO distance, and Hubble constant mea- 
surements, which were used to fix the background, we 



Given the expansion history, the PPF framework |14l . 
[Tsj is defined by three functions and one parameter. 
From these quantities, the dynamics are determined by 
conservation of energy and momentum and the Bianchi 
identities. The defining quantities are g{a,k), which 
quantifies the effective anisotropic stress of the modifica- 
tions and distinguishes the two gravitational potentials, 
f({a), which describes the relationship between the mat- 
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Parameters 


ACDM 


ACDM (with gISW) 


ACDM (with Eg) 




2.235 ± 0.053 


2.232 


2.235 ± 0.054 


2.237 


2.235 ± 0.054 


2.228 




0.1128 ± 0.0036 


0.1137 


0.1124 ±0.0036 


0.1125 


0.1127 ±0.0037 


0.1131 


e 


1.0407 ± 0.0026 


1.0402 


1.0406 ± 0.0027 


1.0406 


1.0406 ± 0.0027 


1.0409 


T 


0.082 ±0.015 


0.079 


0.084 ±0.016 


0.085 


0.082 ±0.016 


0.081 


Us 


0.957 ±0.012 


0.955 


0.958 ±0.012 


0.956 


0.957 ±0.012 


0.957 


HlO^°As] 


3.202 ± 0.038 


3.204 


3.202 ± 0.038 


3.210 


3.202 ± 0.039 


3.202 




0.273 ±0.016 


0.279 


0.271 ±0.016 


0.271 


0.272 ±0.016 


0.273 


Ho 


70.4 ± 1.4 


69.9 


70.6 ± 1.3 


70.5 


70.5 ±1.3 


70.3 


-21nL 


3027.812 


3061.630 


3027.826 



TABLE I: Means, standard deviations (left subdivision of columns) and best-fit values (right subdivision of columns) with 
likelihood for the flat ACDM model using data fi-om WMAP, ACBAR, CBI, VSA, Union, BAG, and SHOES (left column). 
Including one the gISW (middle column) and Eg data sets (right column). 



Parameters 


ACDM (with CA) 


ACDM (with £g&CA) 


ACDM (aU) 




2.228 ±0.051 


2.233 


2.229 ±0.053 


2.229 


2.231 ±0.053 


2.239 




0.1107 ±0.0019 


0.1112 


0.1107 ±0.0020 


0.1112 


0.1107 ±0.020 


0.1108 


e 


1.0401 ± 0.0025 


1.0400 


1.0402 ± 0.0025 


1.0406 


1.0403 ± 0.0026 


1.0408 


T 


0.081 ±0.015 


0.080 


0.081 ±0.015 


0.078 


0.082 ±0.015 


0.078 




0.956 ±0.012 


0.957 


0.956 ±0.012 


0.956 


0.956 ±0.012 


0.957 




3.195 ±0.035 


3.191 


3.195 ±0.036 


3.193 


3.196 ±0.036 


3.185 




0.2642 ± 0.0099 


0.2664 


0.2638 ± 0.0098 


0.2654 


0.2634 ± 0.0098 


0.2622 


Ho 


71.0 ± 1.1 


70.8 


71.0 ± 1.1 


70.9 


71.1 ± 1.1 


71.3 


-21nL 


3032.706 


3032.746 


3066.240 



TABLE IL Same as Table [H but including the CA (left column), both Eg and CA (middle column), and all (right column) 
additional data sets. 



ter and the metric on superhorizon scales, and faio), 
which defines it in the linearized Newtonian regime. The 
additional parameter cr is the transition scale between 
the superhorizon and Newtonian behaviors. For f{R) 
gravity the PPF expressions were developed in [T^] . For 
completeness, we shall review it here. 
At superhorizon scales, the metric ratio 



ffSH(lna) = 



$ ± * 



is obtained from solving the differential equation 
H" B' 



1 



H' 1- B 
ll ^ IP" 



B' 



1 - B 



where k/aH — > 0, and using the relation 



1-B 



(Al) 



$ = 0, (A2) 



(A3) 



where k/aH — > This follows from conservation 

of curvature fluctuation (C = 0) and momentum, con- 
sidering the superhorizon anisotropy relation $ ± 4* = 
BH'q [6|, where q is the momentum fluctuation. 



In the quasistatic regime, we have g^s 
intermediate scales 



g{a, k) 



gsii + 9Qs{cgk/aH)"!> 
1 ± (cgfc/aiJ)"9 



-1/3 and at 



(A4) 



where 



0.71 



B and rig = 2 



IJ. Further, /<; 



-gsH/3, /g = and cr = 1 [ll 

We supply the PPF modified CAMB code [l^ with the 
functions defined above. This gives a good approxima- 
tion for Bq ^ 1. Above this value, the approximation 
begins to break down at intermediate scales and low red- 
shifts, as we have tested by comparing the metric ratio 
predicted by the PPF functions and the exact numerical 
solution for the ordinary differential equations describing 
scalar linear perturbation theory in f{R) gravity for a 
matter-only universe. The expressions defining the f{R) 
scalar linear perturbation theory can be found in @ . This 
deviation shows up in the low multipoles of the CMB for 
high Bq and partially manifests itself when comparing 
constraints on Bq derived from the CMB with the re- 
sults of However, when including gISW measures or 
cluster abundance, the viable values of Bq lie well within 
the regime where the approximation holds. 
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Parameters 


fiR) 


f{R) (with gISW) 


f{R) (with Eg) 




2.223 ± 0.053 


2.206 


2.225 ±0.054 


2.253 


2.224 ± 0.054 


2.206 




0.1123 ± 0.0036 


0.1109 


0.1117 ±0.0036 


0.1133 


0.1125 ±0.0036 


0.1131 


e 


1.0403 ± 0.0027 


1.0392 


1.0403 ± 0.0027 


1.0416 


1.0403 ± 0.0027 


1.0394 


T 


0.083 ±0.016 


0.082 


0.084 ±0.016 


0.090 


0.083 ±0.016 


0.083 


Us 


0.954 ±0.012 


0.950 


0.954 ±0.012 


0.965 


0.954 ±0.013 


0.952 


ln[10^''yl,] 


3.212 ± 0.040 


3.215 


3.209 ± 0.039 


3.200 


3.213 ±0.039 


3.221 


lOOBo 


< 315 


28 


< 43.2 


0.0 


< 319 


30 




0.272 ±0.016 


0.268 


0.269 ±0.016 


0.272 


0.273 ±0.016 


0.279 


Ho 


70.4 ± 1.4 


70.4 


70.7 ± 1.3 


70.7 


70.3 ±1.3 


69.6 


lO'l/flol 


< 350 


46 


< 69.4 


0.0 


< 353 


51 


-2AlnL 


-1.104 


1.506 


-0.696 



TABLE III: Same as Table HI but for f{R) gravity. — 2AlnL is quoted with respect to the corresponding maximum hkeUhood 
flat ACDM model. Limits on Bo and |/flo| indicate the one-sided ID marginalized upper 95% C.L. Note that as Bo ^ 
reproduces ACDM predictions, the slightly poorer fits of f{R) gravity should be attributed to sampling error in the MCMC 
runs. 



Parameters 


f{R) (with CA) 


f{R) (with Sg&CA) 


f{R) (all) 


lOOnb/i^ 


2.209 ± 0.054 


2.204 


2.213 ±0.054 


2.235 


2.216 ±0.054 


2.210 




0.1064 ± 0.0032 


0.1112 


0.1073 ±0.0029 


0.1108 


0.1076 ± 0.0028 


0.1104 


e 


1.0390 ± 0.0027 


1.0398 


1.0392 ± 0.0027 


1.0413 


1.0394 ± 0.0027 


1.0398 


T 


0.077 ±0.016 


0.080 


0.077 ±0.015 


0.084 


0.079 ± 0.015 


0.075 


ris 


0.953 ±0.012 


0.951 


0.954 ±0.012 


0.956 


0.954 ±0.012 


0.951 


ln[10"M,] 


3.175 ±0.0038 


3.209 


3.179 ±0.037 


3.203 


3.182 ±0.0037 


3.193 


lOOBo 


< 0.333 


0.000 


< 0.152 


0.000 


< 0.112 


0.001 




0.247 ±0.014 


0.268 


0.251 ±0.012 


0.261 


0.252 ±0.012 


0.264 


Ho 


72.2 ± 1.4 


70.5 


71.9 ±1.3 


71.4 


71.9 ± 1.2 


70.8 


lO^I/flol 


< 0.484 


0.001 


< 0.263 


0.000 


< 0.194 


0.002 


-2AlnL 


0.802 


0.264 


0.926 



TABLE IV: Same as Table HIl but for f{R) gravity. See also Table Hill 



Appendix B: The iswwll code for f{R) gravity 

We use the publicly available iswwll code [13, [Hi for 
our analysis. Note that we have turned off weak lens- 
ing contributions in the code, focusing only on the gISW 
constraints. The 42 data points of gISW cross correla- 
tions that are used in the likelihood analysis are collected 
from the Two Micron AU Sky Survey (2MASS) extended 
source catalog (XSC) [13, [U, the luminous red galaxies 
(LRGs) and photometric quasars (QSO) of the SDSS [s^ . 
and the National Radio Astronomy Observator y (N RAO) 
Very Large Array (VLA) Sky Survey (NVSS) [5|. They 
are divided into nine galaxy sample bins j (2MASS0-3, 
LRGO-1, QSOO-1, and NVSS) based on flux (2MASS) or 
redshift (LRG and QSO). These data points are a selec- 
tion of multipole bins from all samples, where the selec- 
tion is based on the avoidance of nonlincarities and sys- 
tematic effects from dust extinction, galaxy foregrounds, 
the thermal Sunyaev-Zel'dovich effect, and point source 
contamination to affect the gISW cross correlations [Ol ■ 

In the remainder of Appendix[Bl we discuss the validity 
of the Limber approximation and elucidate the function 



fj{z) that carries information about the redshift distri- 
bution and bias in the context of f{R) gravity. 



1. Limber approximation 

For f{R) gravity the gISW cross correlations are well 
described within the approximations given in Wl A 31 
The applicability of approximating the product DdG/dz 
in Eq. 1^1^ through Eq. ([T5|) and Eq. ([13]) is a direct con- 
sequence of applying the Limber approximation. This 
can easily be seen from Fig. |3] considering the substitu- 
tion k ^ {£ + l/2)/x(z) at the relevant redshifts. The 
accuracy of the Limber approximation itself, in the case 
of ACDM, is at the order of 10% at £ = 2 and drops 
approximately as at higher £ (see, e.g., p^. Is^. Issf) . 
Apart from the multipole the error depends also on 
the width of the redshift distribution, which changes only 
little with f{R) effects. Given the large errors of the cur- 
rently available data points at low £, we conclude that 
the Limber approximation is applicable and furthermore 
very useful since it is numerically faster than an exact 
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integration. 



by the functions 



2. Redshift distribution and bias 

A further ingredient in the iswwll code is the deter- 
mination of the function fj{z). In the Markov chain, 
fj{z) is recomputed when changing the cosmological pa- 
rameter values. The methods by which this function is 
determined differs for each sample, but they are all based 
on galaxy clustering data. 

The 2MASS galaxies are matched with SDSS galax- 
ies in order to identify their redshifts. To obtain the 
nonlinear power spectrum, the Q model for nonlineari- 
ties d^l is applied. Then, the code computes the galaxy 
power spectrum and fits it to measurements, thereby de- 
termining the bias b{z) and Q. Since the required ac- 
curacy for the estimation of bias is only at the order of 
a few tens of percent [l3, this processing is also appli- 
cable to f{R) gravity. The Q model is also adopted for 
LRG galaxies, where the redshift probability distribution 
is inferred with methods described in Ref. [56|. For QSO, 
first, a preliminary estimate for the redshift distribution 
is deduced by locating a region of sky with high spec- 
troscopic completeness, but simultaneously maintaining 
a large area. Taking into account magnification bias 
and fitting bj{z)Ilj{z) using the quasar power spectrum 
and quasar-LRG cross power yields the desired shape of 
fj{z). Note that in f{R) gravity the relationship be- 
tween the metric combination sensitive to gravitational 
redshifts and lensing $_ and the density perturbations 
is modified such that the expression of the lensing win- 
dow function for magnification effects given in Ref. [l7| 
is rescaledby {l + fB.)-^ < (l + /i?o)"' < 1.1 for Bq < 1, 
where (1-1- /a)^^ > 1. Since this is a small effect and 
magnification bias is subdominant to bias and redshift 
distribution, it is save to neglect it. Due to the rather 
low accuracy requirement, we can furthermore adopt the 
ACDM growth factor in the determination of bi{z)Ili{z) . 
Finally, the effective redshift distribution of NVSS is ob- 
tained from cross correlating with the other samples and 
fj{z) is fitted with a F distribution. 



Appendix C: Cluster abundance in f{R) gravity 

The likelihood code of [1^ utilizes the halo mass func- 
tion of Tinker et al. [53 at z = 0.18. which is described 



dn _ pm dhva ^ 



where the distribution /(cr) is given by 



and 



+ 1 



= I P{k)W{kR)k^dk. 



(CI) 

(C2) 
(C3) 



Here, W is the Fourier transform of the real-space top- 
hat window function of radius R. The free parameters 
a, b, and c at z = 0.18 are fitted to ACDM simulations 
in [53 • We define M as the lensing mass, being the more 
conservative a ppr oach than taking it to be the dynamical 
mass (see [Tsl. l58l|). In the large field limit {Bq > 10^^), 
where no chameleon mechanism takes place, these rela- 
tions provide a reasonable fit to f{R) gravity simulations 
as well, using the same values for the parameters as in 
ACDM, and correspond to using the Sheth-Tormen [H^ 
halo mass function wit h sp herical collapse predicted by 
standard gravity (see [l3l l60j). Another approach is 
to alter the spherical collapse calculations, taking care 
of modified forces [H, H^l- Standard spherical collapse 
tends to overestimate the f{R) effects, while the oppo- 
site occurs for the modified spherical collapse calcula- 
tions. Given the constraints that can be achieved from 
the data, we conclude that it is save to assume the large 
field limit, neglect the chameleon mechanism, and use 
standard spherical collapse, i.e., using the ACDM values 
for the free parameters in the Tinker et al. halo mass 
function. While choosing modified spherical collapse is 
the more conservative approach, the standard spherical 
collapse parameters yield a better fit to the simulations 
(see [i^|6^). We remind the reader that the 68% C.L. in- 
ferred from cluster abundance should be taken with cau- 
tion since it approaches the small field limit {Bq < 10~^) 
and intermediate regime, where chameleon effects appear 
that cannot be described through Eq. (jCip . A better fit 
for the halo mass function in f{R) gravity is an objective 
to future work. 
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